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Abstract 

In this work, the benefits of the phase fitting technique are embedded in 
high order discrete Lagrangian integrators. The proposed methodology cre- 
ates integrators with zero phase lag in a test Lagrangian in a similar way 
used in phase fitted numerical methods for ordinary differential equations. 
Moreover, an efficient method for frequency evaluation is proposed based 
on the eccentricities of the moving objects. The results show that the new 
method dramatically improves the accuracy and total energy behaviour in 
Hamiltonian systems. Numerical tests for the 2-body problem with ultra 
high eccentricity up to 0.99 for 10 6 periods and to the Henon-Heiles Hamil- 
tonian system with chaotic behaviour, show the efficiency of the proposed 
approach. 
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1. Introduction 

In the field of numerical integration, methods specially tuned on oscillat- 
ing functions, are of great practical importance. Such methods are needed 
in various branches of natural sciences, particularly in physics, since a lot of 
physical phenomena exhibit a pronounced oscillatory behaviour. For a review 
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of such methods s ee llxaru et al.l (119971); IVandenBerghe et al.l (119991 . l200ll ); 



Ixaru et al.l (120031 ); I VanDaele and VandenBerghd (120070 and references there 



in as well as the book Ixaru and VandenBerg h T(l2004h . 

For problems having highly oscillatory solutions, standard methods with 
unspecialized use can require a huge number of steps to track the oscilla- 
tions. One way to obtain a more efficient integration process is to construct 
numerical methods with an increased algebraic order, although the simple 
implementation of high algebraic order methods may cause severa l problems 
(for example, the existence of parasitic solutions iQuinlanl ( 119991 )). On the 
other hand, there are some special techniques for optimizing numerical meth- 
ods. Trigonometrical fitting and phase-fitting are some of them, producing 
methods with variable coefficients, which depend on v = uh, where oj is the 
dominant frequency of the problem and h is the step length of integration. 
This technique is kno wn as exponent i al (or trigono metric if /i = iu) fitting 
and has a long history iGautschil ( 19611 ). iLychd (1972). An important property 
of exponential fitted algorithms is that they tend to the classical ones when 
the involved frequencies tend to zero, a fact which allows to say that exponen- 
tial fitting represents a natural extension of the classical polynomial fitting. 
The examination of the con vergence of e xponential fitted multistep methods 
is includ ed in Lyches theory iLvche (119721) . The general theory is presented in 
detail in llxaru and VandenBerghd (120041 ) . Furthermore, considering the ac- 
curacy of a method when solving oscillatory problems, it is more appropriate 
to work with the phase-lag, rather than its usual primary lo cal truncation er- 
ror. W e mention the pioneering paper of Brusa and Nigro iBrusa and Nigro 
(119801 ) . in which the phase-lag property was introduced. This is actually 
another type of a truncation error, i.e. the angle between the analytical 
solution and the numerical solution. A significant application of the phase 
or exponential fitting is on the construction of symplecti c methods for os- 



cillatory pro blems encountered in physics and chemistry (jMonovasilis et al. 
(120051 . 120061 )). 

Although phase fitting and exponential fitting are a major improvement 
over algebraic fitted methods especially for oscillatory and orbital problems, 
there is not significant evidence from published results that these methods 
can be applied for long term integration (for example for millions or billions 
of periods). Moreover, several authors use to test their methods to the well 
known 2-body problem but only for relatively low eccentricities (up to 0.2) 
and for relatively small numb er of periods ( no more than several thou s ands). 
We mention here the efforts of lVandeVyverl ( 20061 . 120051 ) ; IWangi ( 20051 ); ISimosI 
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(120041 ) ; lAnastassi and Simoa (120051 . |2004j) in which there is no evidence that 
the phase fitting or trigonometric fitting can be applied to high eccentricities 
(for example to the Halley comet with eccentricity close to 0.967) and for 
long time. 

Another approach to oscillatory and especially Hamiltonian systems is 
the theory of disc r ete va r iational mecha n ics, wh i ch wa s set up in the 1960s 
Jordan and Polakl (ll964T ): ICadz^wl (Il970f ): iLoganl (119731 ) and then it was pro- 



posed in the optimal control literature. It then motivated a lot of authors 
and soon the discrete Euler-Lagrange equations were formulated and the first 
integrators in the discrete calculus of variation and further the multi-freedom 
and higher-order problems were studied. Afterwards, the canonical structure 
and symmetries for discrete system s were o btain e d, and Noether's theorem 



to the discrete case was extended iMaedal (ll980l. 11981). Finally, the time 



as a discrete dynamical variable was regarded iLed (119831 ). A detailed de- 
sc ription of the essential pr o perties of variat i onal i ntegr a tors c an be found 



in 



Marsden and Westl (120011 ); iMarsden et all (119981 ): ILed (119831 ). One of the 



most important properties of variational integrators is that since the discrete 
Lagrangian is an approximation of a continuous Lagrangian function, the 
obtained numerical integrator inherits some of the geometric properties of 
the continuous Lagrangian (such as symplecticity, momentum preservation). 

In the present work, the benefits of the two approach are combined in 
order to construct high order discrete Lagrangian integrators with phase 
fitting. To obtain this, we have adopted a test Lagrangian problem (similar 
to the test ODE in the phase fitting) which is the harmonic oscillator with 
given frequency us. Then, we construct discrete variational schemes that solve 
exactly the test Lagrangian. The application of the method to a general 
Lagrangian needs the determination of the frequency u at every step of the 
integration. The method is applied to the 2-body problem with eccentricity 
up to 0.99 for 10 6 periods and to Henon-Heiles system which for high energies 
exhibit chaotic behaviour. The results clearly demonstrate the efficiency of 
the new approach. 



2. Discrete Variational Mechanics 

The well known least action principle of the continuous Lagrange - Hamil- 
ton Dynamics can be used as a guiding principle to derive discrete integra- 
tors. Following the steps of the derivation of Euler-Lagrange equations in 
the continuous time Lagrangian dynamics, one can derive the discrete time 
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Euler-Lagrange equations. For this purpose, one considers positions g and 
qi and a time step h G R, in order to replace the parameters of position q 
and velocity q in the continuous time Lagrangian L(q,q,t). Then, by con- 
sidering the variable ft as a very small (positive) number, the positions go 
and qi could be thought of as being two points on a curve (trajectory of the 
mechanical system) at time h apart. Under these assumptions, the following 
approximations hold: 

go ~ q(0) , qi ~ q(h) , 

and a function Ld(qo, qi, h) could be defined known as a discrete Lagrangian 
function. 

Many authors assume such functions to approximate the action integral 
along the curve segment between q and g 1; i.e. 

L d (q Q ,q 1 ,h)= L(q,q,t)dt (1) 
Jo 

Furthermore, one may consider the very simple approxi mation for this inte 



gral g iven on the basis of the rectangle rule described in iMarsden and West 



( 120011 ). According to this rule, the integral L Ldt could be approximated 
by the product of the time-interval h times the value of the integrand L ob- 
tained with the velocity q replaced by the approximation (qi — qo)/h: The 
next step is to consider a discrete curve defined by the set of points {qk}k=oi 
and calculate the discrete action along this sequence by summing the discrete 
Lagrangian of the form Ld(qk, qk+i, h) defined for each adjacent pair of points 

{Qki Qk+l)- 

Following the case of the continuous dynamics, we compute variations 
of this action sum with the boundary points go an d g^v held fixed. Briefly, 
discretization of the action functional leads to the concept of an action sum 

n-1 

S d (id) = ^2 L d(ik-i, Qk), id = (go, <?n-i) e Q n (2) 

k=l 

where Ld : Q x Q — > R is an approximation of L called the discrete La- 
grangian. Hence, in the discrete setting the correspondence to the velocity 
phase space TQ is Q x Q. An intuitive motivation for this is that two points 
close to each other correspond approximately to the same information as one 
point and a velocity vector. The discrete Hamilton's principle states that if 
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7(2 is a motion of the discrete mechanical system then it extremizes the ac- 
tion sum, i. e., SSj = 0. By differentiation and rearranging of the terms and 
having in mind that both go and are fixed, the discrete Euler-Lagrange 
(DEL) equation is obtained: 



D 2 L d (q k ~i, q k , h) + D 1 L d (q k , q k+1 , h) = 



(3) 



where the notation DiLd indicates the slot derivative with respect to the 
argument of L^. 

We can define now the map &:QxQ^QxQ, where Q is the space of 
generalized positions q, by which 



D X L A o $ + D 2 L 







(4) 



which means that &(q k -i,q k ) = (q k ,q k +i)- Then, if for each q G Q, the map 
D x L d {q, q) : T q Q -> T*Q is invertible, then D x L d : Q x Q -> T*Q is locally 
invertible and so the discret e flow defined by t he map $ is well defined for 
small enough time steps (see lKane et al.l (119991 ) for details). Moreover, if we 
define the fiber derivative 



FLd-.QxQ^ T*Q 



(5) 



and the two- form u on Q x Q by pulling back the canonical two- form VLqan 
dq i A dpi from T*Q to Q x Q: 



UJ 



FL d {yt C AN) 



The coordinate expression for uj is 



U! 



d 2 L d 
dq{dq{ +1 



Qk+l 



fc+1 



(6) 



(7) 



and can be easily proved that t he map $ preserves the sy mplectic form uj (two 

differe nt proofs are presented in lMarsden et al.l (119981 ) and lWendlandt and Marsden 
( 119971 ) ) . Finally, assuming that the discrete Lagrangian is invariant under 
the action of a Lie group G on Q and £ 6 g, the Lie algebra of G, by anal- 
ogy with the continuous case, we can define the discrete momentum map 
Jd : Q x Q -> g* by 



(Jd(qk,qk+i),0 ■■= (D a L d (q k ,q k+ i,£ Q (q k )) 



(8) 
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It can be proved that the map $ preserves the momentum map J jWendlandt and Marsden 



(11997ft . 

In a position-momentum form the discrete Euler-Lagrange equations ([3]) 
can be defined by the equations below 

p k = -D 1 L d (q k ,q k+1 ,h) 
p k+1 = D 2 L d (q k , q k+1 , h) (9) 

3. Phase-fitted Discrete Lagrangian Integrators 

Summarizing the phase fitting technique, we consider for simplicity only 
first order differential equations, although the same results can be easily 
obtained for second order equations too. Consider the test problem 

^- = iu y(t), y(0) = l (10) 

with exact solution 

y(t) = (11) 

where ujq is a non- negative real value. Let &(h) be a numerical map which 
when it is applied to a set of known past values, it produces a numerical 
estimation of y(t + h). If we assume that all past values are known exactly, 
then the numerical estimation y(t + h) of y(t + h) will be 

y (t + h) = a{u h) ■ e^ot+H^oh)) ( 12 ) 

while the exact solution is e l (^o*+^o^)_ Then the ratio of the estimated to the 
exact solution is 

L = y^lV = a (u h)e-^ h ^ h » (13) 
y(t + h) V ; K J 

In the above equation (fTBl) . the term a(uj h) is called the amplification error, 
while the term 1{uqK) = oj$h — (p(uJoh) is called the phase lag of the numerical 
map. In the case that a(uoh) = 1 and 1{uqK) = 0, we say that the numerical 
map $(/i) is exponentially fitted at the frequency a>o and at the step size h. 
The technique of phase fitting can now be considered as the vanishing or 
minimization of the phase lag. 

Consider now the discrete Lagrangian L d (q k ,q k+1 ,h) (q k corresponds to 
time t k and q k+ i to time t k+ \ — t k + h) and a set of s intermediate points 
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g J with qi = q (t k + cPh). The role of the number of intermidiate points will 
be discussed later. Assuming that c 1 = and c s = 1 we always have q 1 = q k 
and q s = q k +\. Then we can approximate Ld with the quadrature 

s 

L d (q k ,q k+1 ,ti) = h^Wj- L (q(t k + c?h),q{t k + c j h),c j h) (14) 
i=i 

For maximal algebraic order it is easily proved that the following conditions 
must hold: 

xx^y = 7 ^ I , 1=0,1,.. as) 

Then, we can approximate intermediate points and their derivatives with 

q j = b j q k + b j q k+1 , , 

Consider now the test Lagrangian (harmonic oscillator) similar to the test 
equation (TTOT) 

U = \e - (17) 
Then, applying the above assumptions in Eq. ([3]) we get 

qk+1 ~ EUw, (uW - BW^) Qk ~ ^ (18) 
where u = uh. Since the exact solution of Eq. (TTTj) is 

g(t) = Ae iuJt + Se - ^ (19) 

and we want to force our method to solve exactly Eq. ([3]), we get 

IP = cos (c?u) sin (cPu) 

smu 

y sin (c^u) 



smu 

cosu 



B 3 = —usin [c? u) — u— — cos (cPu) 



smu 



cos (c J u 



B 3 = u i - (20) 

smu 
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The role of intermediate points will be clarified now by introducing cor- 
rections at the these points. Let Xj,j = 1, . . . ,S a set of parameters and 
rewrite eq. (flEj) as 



Q 



\>>q k + b> (x s - x S -i) + p(x u x 2 ,..., x S -i, t) 



<? J = \ ( Bj( lk + (x s - xs-i)) + P{xi,x 2 , x S -i, t) 



(21) 



where the p(xi,x 2 , ■ ■ ■ ,xs-i,t) is the interpolating polynomial of the set of 
points {(t k , 0), (t k + c 2 h, xx), (t k + c 3 h, x 2 ), . . . , (t k +i, ^5-1)} and p its time 
derivative. Thus, xi,x 2 , . . . , xs-i can be considered as corrections to points 
q 2 ,q 3 , . . . ,q s (the correction at q 1 is zero since we have assumed that q 1 
The set of equations (jUJ) are now rewritten as 

dL d (q k ,xi,x 2 ,...,x s ,h) 



Qh) 



dq k 

dL d (q k ,xt,x 2 , . . . ,x s ,h) 



+ X(x s - qk+i] 







dx; 



Pk+i 



= 0, j = l,2,...,S-l 
dL d (q k ,x 1 ,x 2 , ...,x s ,h) 

dx.q 



(22) 



where A is the Lagrange multiplier and can be easily prov ed tha t it is equal 
to p k . This technique is similar to those described in Leokl ( 2005 ) and 



Kharevvch et all (l2006h . 



4. Frequency Evaluation 

In order to efficiently evaluate the frequency of the problem, we focus on 
orbital problems and especially on the eccentricity. In general, an elliptic 
orbit may be parameterized as 

X = X m ■ COSU + C\ 

y = y m - sinu + c 2 

where r = (x, y) is the position at time t and u is a function of time. Then 

(23) 

where the product x m ■ y m is equal to the product of the square of the semi- 
major axis a, multiplied by a/1 — e 2 where e is the eccentricity. Since the 



r x r 



Xm ' ym 



du 
~dt 
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frequency locally can be approximated by | ^ | we get 



CO 



r x r 



a 2 VT 



(24) 



Notice here that both a a nd e can b e calcu lated by the position and velocity 
of the moving object (see Goldstein ( 19531 )). 



5. Numerical Tests 

5.1. The 2-body problem 

We now turn to the study of two objects interacting through a central 
force. The most famous example of this type, is the Kepler problem (also 
called the two-body problem) that describes the motion of two bodies which 
attract each other. In the solar system the gravitational interaction between 
two bodies leads to the elliptic orbits of planets and the hyperbolic orbits of 
comets. 

If we choose one of the bodies as the center of our coordinate system, the 
motion will stay in a plane. Denoting the position of the second body by 
q = (qi, q2) T , the Lagrangian of the system takes the form (assuming masses 
and gravitational constant equal to 1) 

L(q,q,t) = Vq + ^ (25) 

The initial conditions are taken 

q=(l-6,0f, q= (26) 

where e is the eccentricity of the orbit. In order to check the efficiency 
of the proposed algorithm, we shall consider only high eccentricities (e = 
0.95 — 0.99). In the first test, we count the number of integration steps needed 
for one period and for eccentricity equal to 0.95. The results are summarized 
in Table [TJ The results have been obtained by adaptively calculating the 
time step, in order to keep the relative error in energy less than 10~ 6 (as 
relative error we mean the absolute error divided by the correct value). 

In the second test, we check the performance of the new method for long 
term integration. First, we integrate the 2-body problem for 10 6 periods 
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Table 1: Number of integration steps 



s 


Linear 


Phase fitted 


1 


11067 


1789 


2 


9873 


1124 


3 


6534 


252 


4 


4321 


145 


5 


3245 


46 



and for eccentricity 0.99. Figure Q] shows the exact orbit (solid line), the 
calculated points for the first period (o) and the calculated points for the 
last period (□). Again the time step is adaptively calculated in order to keep 
the relative error in energy less than 10~ 6 (the value 10 -6 of course can be 
changed to obtain higher or less accuracy, but in these tests has been selected 
because it produces 10 full periods in less than a second in a typical personal 
computer). Figure [2] shows the solution produced for the perturbed Kepler 
problem described by the Lagrangian 

L(q,q) = -q T q+ — + -—^ (27) 
2 |q| 2|q| c> 

again for 10 6 periods and for eccentricity 0.6 where it is clear that the solution 
is an ellipse that rotates slowly around one of its foci. Again the time step 
is adaptively control in order to keep the relative error in energy less than 
10~ 6 . All the previous test use S = 5 as the value of intermediate points. 

5.2. Henon-Heiles Hamiltonian System 

In second test, we examine the behaviour of the new method in the Henon- 
Heiles Hamiltonian system. In the 1960s, a model of the motion of stars in 
a cylindricall y symm e tric, t ime-independent pote ntial were investigated b y 
astronomers ( Vernov ( 2003 )). Henon and Heiles ( Henon and Heiles ( 19641 )) 



proposed the Hamiltonian 

H = \ {x 2 + y 2 ) + \ [x 2 + y 2 ) + x 2 y - ^y 3 (28) 

where for small values of energy the trajectories are trivial but for higher 
energies, dynamic chaos emerges in the system. Setting the total energy as 
E = 2c 2 , Petrov produced asymptotic solutions of the form of the product of 
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Figure 1: The exact orbit for the 2-body problem for eccentricity 0.99 for 10 6 period (solid 
line), the calculated points for the first period (o) and the calculated points for the last 
period (□). 



1.5 




-f.5 ^1 ^5 05 1 1.5 

X 



Figure 2: The calculated orbit for the perturbed Kepler problem for eccentricity 0.6 and 
for 10 6 periods. 



11 



one s low and one fast oscillation with a characteristic period T (see iPetrov 
(120071 )). In our test, we calculate the winding number of the orbit around 
the origin for a half-period. The calculated values are compa r ed wi th the 
theoretical ones produced by the asymptotic solution (Petroy ( 20071 )) and 
with a set of very accurate values produced in the following way: first the 
RKN86 8-stages Runge-K u tta-N ystrom pair of orders 8 and 6 was used (see 
Papakostas and Tsitouras (Il999f a The error tolerance of the method was 
set to IE — 14 almost at the machine precision. Then, at each step, the 
calculated solution was projected in the manifold described by the equations 
E = const and J = const where E, J the total energy and angular momen- 
tum respectively. Again the method that was applied uses S = 5 and an 
adaptive time step calculation keeping the relative error in energy less than 
10~ 9 . Table [2] presents the results for c = 0.1, 0.05 and 0.0025. Notice here 
that for c = 0.1, dynamic chaos is present in the system. 



Tabic 2: Winding number for half-period for the Henon-Heiles System 



c 


Theoretical 


Our Method 


RKN86 method 


0.1 


430 


425 


425 


0.05 


1700 


1697 


1697 


0.025 


6800 


6865 


6865 



6. Conclusions 

It has been shown in this work, that the technique of phase fitting, when 
it is embedded in discrete Lagrangian integrators, improves the accuracy and 
the efficiency of the numerical method. Following the classical application 
of the phase fitting technique, the discrete Lagrangian integrator is forced to 
solve exactly the test Lagrangian of harmonic oscillator with a given self fre- 
quency. The coefficients of the resulting integrator, depend on the frequency 
of the problem at each integration step. In order to improve the accuracy 
of the method, a set of intermediate points were added as corrections to the 
trigonometric path. The results show that the method can be used for long 
term integrations of planetary motions and this was demonstrated by apply- 
ing the method to very high eccentricities (0.99) and for milions of periods. 
Moreover, a simple but quite efficient technique for frequency evaluation is 
proposed based on the eccentricity of the integrated orbit. 
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